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ABSTRACT 

We analyze stellar convection with the aid of 3D hydrodynamic simulations, intro- 
ducing the turbulent cascade into our theoretical analysis. We devise closures of the 
Reynolds-decomposed mean field equations by simple physical modeling of the simu- 
lations (we relate temperature and density fluctuations via coefficients); the procedure 
(CABS, Convection Algorithm Based on Simulations) is terrestrially testable and is 
amenable to systematic improvement. We develop a turbulent kinetic energy equation 
which contains both nonlocal and time dependent terms, and is appropriate if the con- 
vective transit time is shorter than the evolutionary time scale. The interpretation of 
mixing- length theory (MLT) as generally used in astrophysics is incorrect; MLT forces 
the mixing length to be an imposed constant. Direct tests show that the damping as- 
sociated with the flow is that suggested by Kolmogorov (ex ~ pW)^ms/ ^d, where Id 
is the size of the largest eddy and {u')rms is the local rms turbulent velocity). This 
eddy size is approximately the depth of the convection zone £cz in our simulations, and 
corresponds in some respects to the mixing length of MLT. New terms involving the 
local heating due to turbulent dissipation should appear in the stellar evolutionary equa- 
tions, and are not guaranteed to be negligible. The enthalpy flux (stellar "convective 
luminosity") is directly connected to the buoyant acceleration, and hence to the scale 
of convective velocity. MLT tends to systematically underestimate the velocity scale, 
which affects estimates of chromospheric and coronal heating, mass loss, and wave gen- 
eration. Quantitative comparison with a variety of 3D simulations reveals a previously 
unrecognized consistency. Extension of this approach to deal with rotational shear and 
MHD is indicated. Examples of application to stellar evolution will be presented in 
subsequent papers in this series. 

Subject headings: stars: evolution - hydrodynamics - convection - turbulence 



^Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson AZ 85721 
"^FLASH Center, University of Chicago, Chicago, IL 

^School of Earth and Space Exploration, Arizona State University, Tempe, AZ 
*Joint Institute for Nuclear Astrophysics, University of Chicago, Chicago, IL 



- 2 - 



Introduction 



More than fifty years ago, the version of the "mixing length theory" of convecti on (MLT) which 
became the p r eferred basis for subs equent study of stehar evolution was introduced ( Biermannlll951 



Vitense 



1953 



Bohm- Vitensd 1 1 958l ) . Despite much effort (still ongoing), MLT is still the standard 



choice for the field. 

In this paper we develop a new procedure, "Convective Algorithms Based on Simulations" 
(CABS), in which we close the Reynolds-decomposed, angular and time averaged equations by sim- 
ple physical models based upon analy sis of fully three-dimensional, time dependent turbulent stellar 



convection (jMeakin &: Arnettll2007bl ). These simulations include convective boundaries within the 



computational volume (as far from the edges of the grid as feasible), and allow interface physics 



198C 



simula t ions (especiallv IChan &: Sofia 
(|200d ): IPorter. Woodward, fc Jacobsl (l2000l) : 



1996): 


Kim et al. 


(1995. 


1996): 



Robinson et al 



(1200411). and from research in other 



fields , such as terrestrial fluids (lTurnerlll973l ). oceanography (lGilllll982l ). and meteorology (jPutton 
19861 ). which have a firmer empirical basis. We develop a simple description of the convective 
velocity field as seen in our simulations. This effort brings some startling suggestions for revision 
of our interpretation of M LT, and suggests how our approach may be generalize d to include ro- 
tation and magnetic fields (jBalbus k, Hawlevlll998l : iPessah. Chan, k, Psaltisll2006l ) . This is timely , 
considering recent succe ss in simulating turbulent plasma with magnetic fields (|Browningl 120081 : 
Schiissler k Vokeilboosl ). 



Since the formulation of MLT, there has been a considerable devel opment in understanding 
the nature of chaotic b ehavior i n non li near sy s tems; see ICvitanovicI (Il989l) for a rev ie w and r e print s 
of original papers, and iFrischI (|l995l ): ICleickl (jl987l ): [Thompson k Steward! (119861). iLorenzj (jl963l ) 



1961 



which 



presented a solution to the Rayleigh problem of thermal convection (| Chandr asekharl 
captured the seed of chaos in the Lorenz attr actor, and contains a re presentatio r i of th e fluctuat- 
ing aspect of turbulence not present in MLT. iKolmogorovl (|l962l ) and lObukhovl (|l962l ) developed 
the modern versio n of the turbulent cascade. Although already formulated, the original theory 
(jKolmogorovlll94ll ) was not used in MLT. 



We derive our approximate theory from a consideration of the full equations of 3D compressible 
hydrodynamics for a multiple component fluid. These are close to the corresponding equations for 
a high beta (matter dominated) plasma. This approach will allow us to incorporate a variety of 
phenomena in a coherent way, and to evaluate their relative importance, rather than to patch 
together various bits of physics piecemeal. Here we focus on the dominant features of non-rotating, 
non-magnetic, turbulent, compressible fluid flow. The results are applicable to almost all stages of 
stellar evolution (any stage having convection or shear )0. 



^The shear from convection is similar to the shear from differential rotation; fluid experiments may use either to 
investigate the physics of shear ( Turneij|l973l ). Although different in some details, there are deep connections between 
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In Section 2 we examine the physical aspects of stellar convection, and show that the velocity 
scale is set by the balance between buoyant driving, and damping in the Kolmogorov cascade. 
Appropriate averaging allows us to deal with mild time dependence and reveals a robust underly- 
ing behavior. We develop a kinetic energy equation describing the average properties of turbulent 
convection, which shows how turbulent motion is created, transported, and destroyed. In Section 3 
we incorporate this theoretical development into the equations of stellar evolution with turbulent 
convection, including new terms. Section 4 uses the theoretical development to compare our simu- 
lations to others, and finds previously unrecognized similarities. Section 5 indicates some important 
implications of this work, including how the effects of burning, rotation and magnetic fields may 
be included. Section 6 summarizes our results and conclusions. Quantitative treatment of the 
dynamics of fluctuations and a comprehensive algorithm for stellar evolutionary calculations will 
be developed in subsequent papers in this series. 



2. Physical Aspects of Stellar Convection 

Stellar convection has high Reynolds numbers because of the large linear scales, and is therefore 
highly turbulent. Our simulations have adequate resolution to show this type of behavior. Turbulent 
convection has several key features that need to be modeled: (1) it is nonlocal, (2) it has strong 
fluctuations in both space and time, (3) it is damped by a cascade of unstable vortices down to 
scales small enough for microscopic processes to dissipate effectively, (4) mixing of passive scalar 
properties is efficient, (5) turbulent behavior spreads to fill the volume allowed, and (6) buoyant 
acceleration is closely related to convective (enthalpy) flux, so that convective kinetic energy is 
closely tied to convective luminosity. MLT incorporates (4) and imperfectly deals with (6); we will 
address all six issues. 



Cattaneo. et al. I (|l99ll ) found that 3D simulations of turbulent convection had two aspects: 
(1) vigorous, large scale downflow^, and (2) disorganized weaker motions. These two aspects 
have been confirmed by many other simulations, including our own. MLT attempts to describe 
the average properties of the disordered aspect. We will construct a theory which includes both; 
buoyant acceleration is characteristic of the largest scales, and turbulent dissipation of the smallest. 



2.1. The Kinetic Energy Equation 



We start with the equations we used in our 3D simulations (jMeakin &: ArnettI l2007bl ). We 
use Reynolds decomposition of relevant flow and thermodynamic variables , which separates the 
fluctuating from the slowly varying components (iTennekes Lumleyi Il972l ). Numerical simula- 



convective mixing as described here and the rotational mixing investigated by Mevnet fc Maedei 1 2000l ). 
^Porter fc WoodwardI (|200ol ') suggest that large scale flows do dominate the energy flux. 
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tions dissipate features with wavelen gths at or below the grid scale, and we will identify this with 
dissipation in the turbulent cascade (|Kolmogorovl Il94ll . Il962l ) . 



In the process of approximation of partial differential equations by finite methods, there is 
inevitably a loss of information at scales smaller than the grid size. A single volume element in 
space is approximated as a homogeneous entity; this is equivalent to complete mixing at this scale, 
at each time step, of mass, momentum, and energ y. The loss of information that occurs with 



this mixing corresponds to an increa se in entropv (IShannon 



equivalent to the action of viscosity (iLandau &: Lifshitall959l ). In 3D flow, turbulent energy will 



19481 ): the mixing of momentum is 



cascade from large scales to small, at a rate set by the largest scales. We use an implicit sub-grid 
dissipation in our large eddy simulation (ILES), which is the most computationally e f ficieri t way 
to deal with turbulent systems with a large range of scales (jBorisI 120071 : IWoodwardI 120071 ) ; the 
largest scales, which set the rate of cascade and contain most of the energy, are resolved on our 
grid and explicitly calculated, while the sub-grid scales are dissipated in a way consistent with the 
Kolmogorov cascade. 



Svtine. et al. I (|2000l ) demonstrated that PPM, the piece-wise parabolic method based on the 



Euler equation (which has no explicit viscosity), converges to the same limit as methods based on 
compressible viscous equations (which do have explicit viscosity), as th e grid is refined to s maller 
zones and smaller effective viscosity (the relevant lim it for astrophysics). iPorter et al. I (119991^ show 



the co mpatibility of mildly compressible flow with the lKolmogorovl ()194ll ) spectrum: iKritsuk. et al. 



(| 20071 ) have pushed this to highly compressible flows as well. To represent the sub-grid dissipation, 
which is inherent in our simulations, we explicitly introduce a volumetric dissipation rate for kinetic 
energy, ek, in our theoretical analysis. However, the turbulent cascade is a property of the whole 
convective flow, so connection to turbulence theory must be made through integrals of ex over the 
convection zone. N o te that this is differe n t from defining ek a s a function of local variables (e.g., 
Smagorinsk^ (|l963l l: Ichan fc Sofial (|l989l ) : iHansen k Kawaleil (|l994l ^: ksida ArnettI ^20od )): see 
below. 



The kinetic energy (KE) equation for convective motion was given in lMeakin Sz ArnettI (|2007bl ). 
Taking the scalar product of the velocity with the equation of motion, we decompose the convective 
velocity u, the density p, and the pressure p into mean and fluctuating components (e.g., p = Po+p', 
so the time averages are p = Po and p' = 0. This choice of just u, p, and p for this Reynolds 
decomposition into average and fluctuating parts gives the simplest equation for kinetic energy; the 
velocity u is derived from buoyancy and pressure forces {p and p fluctuations). 

Using the hydrostati c equilibrium condition, and performing averages, gives (see eq. A. 12 in 
Meakin k ArnettI (|2007bl ^). 



dt{pEK) + V ■ (pEkmo) 



-V • (Fp + Fk) + {p'V ■ u') 

+ (/9'g • U') - £k. 



(1) 
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We use {p) to denote an average over angles at constant radius. The time average is taken over 
durations greater than a transit time ttransu = '^cz/{u')rms, where icz is the depth of the convective 
region, and {u')rms is the rms velocity across this region. This smooths the fluctuations, gives a 
nonlocal character to the analysis, and implies a separation of time scales into short (t <C ttransu) 
and long {t S> ttransit)- We consider the case in which we may integrate over these short time scales 
and explicitly calculate the evolution on the long time scales. 

Here Ek is the kinetic energy per unit mass, p the mass density, uq is the nonfluctuating part 
of the fluid velocity vector and u' its fluctuating part, so 

Fp = p'u' (2) 

is the energy flux due to pressure perturbations carried by fluctuations (thi s pressure-velocity 



correl ation flux reduces to the acoustic flux when considering sound waves, iLandau &: Lifshitz 



(|l959l ). p. 251, and contains the energy flux due to internal gravity waves), 

Fi^ = pEku' (3) 

is the flux of kinetic energy carried by convective turbulent motion, g is the gravitational accelera- 
tion vector, and £k is the volumetric rate of dissipation implied by the turbulent cascade down to 
small scales. 



2.2. Boundaries 

The location of convective boundaries is a long-standing problem in stellar astrophysics. It 
has long been known that some sort of " convective overshooting" is necessary for models to match 
observations. This has been conceived as mixing of material beyond the convective boundary, often 
by "penetrative convection". This problem exis t s in p art because the convective boundary has 



been inappropriately deflned (iMeakin &: ArnettI (|2007bl ) . especially §4.1 and §7). The standard 



deflnition used in mixing length theory defines the convective boundaries using the thermodynamic 
nablas, which essentially mark the onset and cessation of buoyant acceleration. One of the essential 
properties of convection is that turbulence fills the space available. A more appropriate criterion 
is one in which the stellar background is stiff enough to contain the turbulent convection, and 
therefore must account for the relative strength of the turbulent fiow and the elastic stable layers. 

We will define the convective zone as that region in which the stratification of the medium is 
unstable to turbulent mixing. This is evaluated with the "bulk Richardson number" 

RiB = Ab l/u^, (4) 

where Ab = J^^ N'^dr is the buoyancy jump across a layer of thickness Ar in the radial direction, 
N"^ is the Briint-Vaisala frequency, u is the rms velocity providing shear, and I is the scale length 
of the turbulence (essentially the size of the largest eddies) . The precise definition of the thickness 
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Fig. 1. — Behavior of global quantities in the convection zone, which affect the turbulent kinetic 
energy, entrainment, and gravity wave generation at the convective boundaries. The thick line is the 
sum of all terms exce pt entrainment and boundary wave luminosity, and indicates when entrainment 
events are vigorous (jMeakin fc ArnettI (|2007bl ). Fig. 4) . The buoyancy driving B = J W^dV is 
compensated for by the turbulent damping D = J e^dV. The time derivative of the kinetic energy 
in the convection zone is denoted dEx/dt = ^ / p{u')1^gdV . There is a slow increase in amplitude 
due to lack of balance between nuclear heating and neutrino cooling (see Fig. [5] below and lArnett 
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Ar is a topic deserving further study; in several cases we have noticed that the rapid change in N"^ 
near a boundary tends to make A6 insensitive to the exact value of Ar chosen for integration to 
obtain A5. 

In the linear limit of plane par allel flow, a re gion with Rig < 0.25 has enough kinetic energy to 
overcome the stable stratification (lDuttonlll986l ). Here Rig is the "gradient" Richardson number, 
and is a locally defined quantity, used in a linear stability analysis. The "bulk" Richardson num- 
ber Ris is an inherently nonlinear quantity based on integration over an extended region. Since 
Ria crit > 0, this allows for N"^ > as well. This formulation recognizes that layers that are ther- 
modynamically stable (real roots to A^^) can be hydrodynamically unstable when kinetic energy is 
input to the zone. The bulk Richardson criterion allows more mixing than predicted by the Ledoux 
criterion; the Sch warzschild criterion igno r es co m positional effects, and als o predicts more mixing 
than Ledoux (see iKippenhahn &: Weigerti (jl99Cl ): iHansen k. Kawalerl (j 19941 ). 



Stable regions adjacent to the convective zone will become Richardson unstable periodically 
as shear builds up from adjacent turbulent motions and waves gen erated by convection, leading to 



entrainment of material at convective boundaries (| Fernando! 1 1 9 9 ll ) . The bulk Richardson number 
is not only an indicator of a convective boundary; it also determine s the "rate" at which the 



boundary migrates through the lagrangian mesh by mixing processes (jMeakin fc ArnettI l2007bl ) 
and thus provides a "dynamic" boundary condition for the flow. 



2.3. Averages 

Our analysis is based primarily upon the numerical simulations of oxygen shell burning in a 
23M0 star, but is found to be broadly applicable to other examples of stellar convection (such 
as convective envelopes, which are driven by the superadiabatic layer at the top rather than the 
nucl ear burning shell a t the bo ttom) . Our convective zone has a depth of two pressure scale heights 



fsee iMeakin fc ArnettI (|2007bl ^ for detail). 



Why do we need to do averaging? Fig. [T] shows the behavior of the largest terms in the kinetic 
energy equation (averaged over the convection zone), throughout the duration of the simulation. 
The dominant terms are the integrals of the buoyancy driving J" W^dV = / gp'u'dV, which is 
positive, and of damping, J exdV, which is negative. Both show recurrent bursts (ten in 800 
seconds). The time derivative of kinetic energy in the convection zone (labeled dEx/dt) is constant 
on average. The remaining bumps are due to missing terms (divergence of Fp and Fx)- 

We saved complete data every 0.5 seconds, which was adequate for making movies, and con- 
structing averages of motion. A key issue is distinguishing between turbulent velocities and wave 
velocities with such a coarse stride in time. We have found a way to split the velocity averages into 
turbulent and wave components for better accuracy in the analysis (see Appendix A for a detailed 
discussion of how this is done). 



-8- 



Table 1. Kinetic Energy Equation Terms Averaged over 
Convection Zone and over four transit times (erg/s) 



Term 


Value 


Term 


Value 


J{WE)dV 


4.576(45) 


-J{W)dV 


-4.677(45) 


- J{V-FK)dV 


2.584(44) 


-/(V • Fp)dV 


-9.922(43) 


-J{dEK/dt)dV 


-5.790(43) 


J{p'Vu')dV 


-1.516(41) 
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The level of complexity shown in Fig. [T] is daunting, but the near balance of buoyant driving 
and turbulent damping provides a clue. 

The system is simplified if we integrate over time; the resulting values over four transit times 
are shown in Table [H We do not calculate the damping directly, but deduce it as the remainder left 
from the other five terms, so that the entries sum to zero. We are comfortable with this procedure 
because of the good conservation properties of the numerical simulations. 

The dominance of buoyant driving and damping is now clear; the next largest term is the 
divergence of the kinetic energy flux, at less than 6 percent of the buoyancy term. 

Suppose the globally-averaged damping term has the Kolmogorov form, so that 

exdV = MczVrmsV^D, (5) 

cz 

where is a constant "damping length," and Vrms is the rms velocity determined over the 
entire convection zone. For a turbulent spectrum, is taken to be the largest length scale 



([Landau &: Lifshitd Il959l ). In what follows we will use the term Kolmogorov damping to mean 



damping due to a turbulent cascade to small scales, having this cubic dependence on rms velocity. 

Unlike MLT, we relate the damping length to the largest scale in the flow; this may be thought 
of as a "coherent structure" or a "plume which traverses the depth of the convection zone." It is 
not a free parameter. We introduce the notation 

an = ^d/^cz, (6) 

so that Q£) woul d be constant i f the s ize of the largest eddy scales simply with the depth of the 
convection zone. iKritsuk. et al. I (|2007l ) found the Kolmogorov damping to be very close to constant 
during a statistical steady state in their 3D (2048^) simulations. 



-10- 



Table 2. Time Scales 



Time scale 


Definition 


Value(seconds) 


Buoyant rise 


/ WbdV/KE 


14.1 


Velocity damping 




27.6 


KE damping 


iD/{2Vrms) 


13.8 


Transit time 


(■Cz/Vrms 


51.4 


Turnover time 


'^(■Cz/ Vrms 


102.8 
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The mass contained in the convection zone is Mqz = 1-84 x 10^^ g and the total kinetic energy 
is J ExdV = 8.61 x 10^^ ergs so that the rms velocity is Vrms = 9.66 x 10® cm/s. This gives 
ijj = 3.56 X 10^ cm « 0.85 icz- It is natural to compare the damping length to the depth of 
the convection zone. Not only is icz the largest length available for an eddy, but if measured in 
pressure scale height units {icz /Hp), it indicates the degree of thermodynamic anisotropy across 
the convective region. 

Table [2] gives several relevant time scales in seconds. Although these times are short in human 
terms, they are much closer to thermal relaxation times than are the corresponding numbers for 
deep simulations of the solar convection zone (our simulations are much more relaxed in this sense) . 

We may write the global damping as MczV^^g/ in = \MczVrms/ i^D I'^Vrms)^ and see that the 
damping time for kinetic energy in the convective zone is half the time to transit a damping length 
(which is approximately the depth of the convection zone) . The turnover time for the convection 
zone is 2r, where r = icz/vrms is the transit time, to be precise. The rise time for kinetic energy 
and the corresponding damping time are similar (14 seconds), and much shorter than a turnover 
time. 

The term "convective efficiency" in the stellar context is usually taken to mean that the time 
scale for convective energy transport is much less than the time scale for radiative energy transport 



(jHansen &: Kawaleii (| 19941 ). p. 187); this insures that the actual temperature gradient is only slightly 
in excess of the adiabatic one. Convection can also be thought of as a thermodynamic cycle, taking 
a time 2icz /vrms- As we saw above, the dissipation time scale for kinetic energy is about one- 
seventh of this (0.134), so that in this sense, convection is not thermodynamically efficient at all, 
but requires continual work to keep it running. The two uses of "efficiency" causes confusion, 
and downplays the fact that stellar convection is highly dissipative, even for slightly super adiabatic 
temperature gradients. 



2.4. Anisotropy and Kinetic Energy 

Here we present an initial, qualitative discussion of the structure of the fl ow, which we are 



i nvest igating quantitatively for subsequent publication. We extend the idea of ICattaneo. et al. 



(|l99ll ) that the convective velocity field has two components, a more ordered global flow and a 
chaotic turbulent flow. We focus on the source and sink for convective kinetic energy. Gravita- 
tional acceleration breaks the symmetry of space; we choose our z-axis parallel to this acceleration. 
Buoyant acceleration starts an anisotropic flow in the z direction. The flow is unstable, and begins 
to break up into smaller scales, and becoming more isotropic. 

Suppose that the flow occurs in narrow, vigorous down plumes and slow wide upflows for 
convective zones with significant anisotropy. The kinetic energy in such a downflow would dominate 
that in the upflow; the kinetic energy flux has two more powers of velocity than the mass flux, so 
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Fig. 2. — Comparison of terms in the kinetic energy equation, averaged over four transit times 
and over angle, as a function of radius. The solid line is the inferred local value of the volumetric 
dissipation due to the turbulent cascade. The thick line of triangles represents a term in the equation 
for kinetic energy density (Eq. [1]). Buoyancy driving (left) and ek approximated as pu^^g/lo = 
1.54/9M^/££) (right) are shown. Here we distinguish between the global rms convective velocity VT-^ns 
and the rms convective velocity at a given radius, {u''^)^ = {u')rms- We find ek = (vrms)^ /^D for 
£d = 3.5 X 10^ cm ~ £cz, globally, and with this choice of damping length, the local damping per 
unit mass is proportional to u^^g/iD- Much of the variati on shown in the bo t tom pa nel is simply 



due to the density gradient (see the angular velocities in iMeakin &: ArnettI (|2007bl ). Fig. 6, left 
panel, and the background density structure in their Fig. 2, top left panel). Note that, as shown, 
the signs of the terms are consistent with the kinetic energy equation, so that up means increasing, 
down decreasing kinetic energy for the terms indicated by triangles. 
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that the fastest flows dominate the kinetic energy budget ([Meakin &: Arnettll2007bl ). We separate 
the kinetic energy into two components, Ma/2, which corresponds to the largest eddies and 
which represents all of the turbulent cascade. The turbulent instability begins to make the flow 
more isotropic, so that 

ul^ul + ul, (7) 



in the z-direction (see Ch. VI in iBatchelorl (ll96Cll )). If this kinetic energy is shared with the two 
perpendicular directions, then 

(8) 



9 2 



Further, it seems that the cascade component in the z-direction is similar to the components in 



the X- or the y-directions, 



u 



b- 



While the downflow matter is balanced by a return flow 



in X and y, the return flow is slower and broader (by our assumptions), with lower kinetic energy 
(which we neglect here). Thus 

ul « ul (9) 

One quarter of the kinetic energy is in the largest scale component u'^/2. The ratio of vertical to 
horizontal rms velocities is therefore u^/uz ~ l/\/2. 



How does this simple model compare with simulations? In Fig. 6, iMeakin &: ArnettI (|2007bl ) 
give the rms values of velocity in the vertical {vr) and the horizontal {vq and f^) directions. We 
identify locally the r direction with the z-axis, and theta and phi with the x and the y axes. Then 
from their Figure 6, we find 

0.5 < Ux/uz ~ Uy/uz < 0.8, (10) 

over the convection zone, away from the boundaries. The measured value of the vertical velocity 
contains both large scale and cascade components. This is consistent with our estimate of u^/uz ~ 
Uy/uz ~ l/\/2 ~ 0.7 07. The large sca l e coin ponent of kinetic energy is comparable to the cascade 
one, consistent with ICattaneo. et al. I (|l99ll ). and our simulations. 



The simulations have led us to a simple, two-component model for the origin and destruction 
of convective kinetic energy. Anisotropy is due to buoyant acceleration in the largest eddy, which 
represents one component. The second is isotropic turbulence, which gives dissipation of kinetic 
energy after the cascade to small scales. It will be interesting to test and refine this model against 
a wide variety of simulations. The power spectra of the flow will be a useful tool to explore this 
ansatz. 

In what follows we will define a turbulent rms velocity ut from the measured horizontal velocity 



(1)2 = 1.540. 



components, u^ = |(m^ + u^) = \u^rns- We note for future reference that m^^^/m; 
Roughly speaking, we identify ut with the turbulent cascade. The anisotropic component Ua is 
identified predominately with the largest turbulent scale, and may be sensitive to the structure of 
the convection zone and the location and strength of the buoyant acceleration. 
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2.5. Radial variation of Averages 

We now relax one level of averaging, and consider the radial variation of all the terms in the 
kinetic energy equation. Figures [2l [3] and U] show six terms in six panels, each compared to the 
inferred damping (the solid line present in each panel). This allows us to determine what each of 
five terms contribute to the inferred damping, and also shows an isotropic Kolmogorov estimate as 
a simple approximation (Fig. [21 right panel) . The solid line has a sharp spike at the bottom of the 
convection zone. This feature is due to the dominance of g-mode character in the motion near the 
boundary region. The Kolmogorov cascade is appropriate away from the boundaries. 



In Figure [21 the left panel shows the averaged buoyancy J {WbcIV), and the right panel shows 
the isotropic Kolmogorov approximation, both as dotted curves. In planar geometry (a fair approx- 
imation), we can relate this q uantity to the "buoyancy flux" which is widely used in experimental 



fluid mechanics (|Turneiill973l ). The buoyancy flux is 



q = (g • u'p')/po, (11) 



so that J (WBdV) f {q)4'Kr'^dr. The enthalpy (convective heat) flux is 



Fe = {poWCpT), (12) 



where Cp is the specific heat at constant pressure. For low Mach number flows like ourqfl the 
pressure fluctuations are small. The temperature and density perturbations at constant pressure 
are proportional, 

p'/po = Pt T'/To, (13) 



where Pt = —dlnp/dlnTj , taken at constant pressure (and composition We find, using 
hydrostatic equilibrium of the unperturbed star Hpg = Pq/po, 

Fe = poHpq/Va. (14) 

where Va = PtP/ pCpT, a factor that we will see again. Except possibly for extreme cases, the 
enthalpy flux Fg, which is intimately related to the stellar lumin osity, is itself close l y relat ed to the 



buoyancy flux q and hence to the convective velocity (see also iMeakin &: ArnettI (|2007bl )). MLT 
ignores this connection; we shall exploit it. The source of turbulent kinetic energy is directly 
proportional to the convective luminosity and thus to the radial entropy gradient ("superadiabatic 
gradient") of MLT. 



■^As we discuss below, simulations with stiffer equations of state will have larger pressure fluctuations. We expect 
the neglected terms to be important for wave generation, but probably in a restricted volume of the convection zone. 

'^In iMeakin fc ArnettI (|2007bl ) we denoted this quantity by /3t; the subscript T is to avoid confusion both with 
Eddington's use of /3 for the ratio of gas to total pressure, and the fluid dynamics community use of /? as the adverse 
temperature gradient. 
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Buoyancy is one of the two dominant terms, and with a sign opposite to that of the damping; 
it differs in magnitude from the inferred damping at the convective boundaries (lower) and the 
broad peak above the bottom of the convection zone (higher). 

The right panel shows the isotropic Kolmogorov approximation to the local damping, Mcz u^msl^-D = 
1.54 Mczu^/(-D- The condition of global balance between all driving and damping terms gives a 
value for but the local value of ut is used. This gives a relatively good, smoothed fit to the actual 
inferred damping throughout the turbulent convection zone, and departs only in the boundary lay- 
ers where the velocity field is due primarily to gravity waves. The net effect of the remaining terms 
in the KE equation is to modify the velocity field so that the damping is more smoothly distributed 
over the turbulent region than is the driving. 

The agreement in FigureEl right panel, between the actual inferred dissipation and our estimate 
provides strong support for our introduction of the "Kolmogorov dissipation" in both its global and 
local forms. 

Fig- El shows the spatial behavior of the flux divergence terms, V • Fk for flux of kinetic energy 
and V • Fp for the pressure correlation flux. These terms have signiflcant positive and negative 
contributions, which cancel upon averaging over radius, and therefore are more important than 
Table [U would suggest. The change of sign in a divergence implies that these terms move kinetic 
energy. In particular, they remove kinetic energy from the region where buoyancy is strong, and 
add it to regions where the buoyancy is weak. The divergence of the kinetic energy flux (Fig. O 
left panel) is the most signiflcant in transporting energy. It moves kinetic energy from the region 
in which buoyancy driving exceeds the inferred damping, and toward the convective boundaries. 

In Fig. [21 left panel, we saw that there is negative buoyancy at the convective boundaries; 
this is due to buoyancy braking of the convective motion. The pressure correlation flux (Fig. [3l 
right panel) is most effective in the deflcit regions right at the convective boundary, and generates 
elastic response (waves) in the stably stratifled regions outside the convection zone. This two-step 
behavior, with flrst Fke and then Fp carrying energy to the edge of the convection zone, is mirrored 
at the upper convective boundary, but at lower amplitude. 

The p'V • u' term, shown in the left panel of Fig. [H does the same sort of thing at a much 
reduced level. In the right panel is shown the time derivative of the kinetic energy, which is small. 
The convection is close to a steady state behavior. 

There is a simple picture which explains these trends. (1) The extent of turbulence is limited 
by energy and by boundaries. In a star, turbulence will mix even stably stratifled layers so long 
as sufflcient kinetic energy is available to supply the work necessary. (2) Turbulence takes ordered 
motion on the large scales and converts it to disordered motion on small scales. It makes the 
flow more isotropic. Kolmogorov dissipation is derived with the assumption of homogeneity and 
isotropy, and so becomes a better approximation as turbulence acts. 

The convective motion is driven by buoyant acceleration, parallel to the gravitational acceler- 
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Fig. 3. — Comparison of terms in the kinetic energy equation, averaged over four transit times 
and over angle, as a function of radius. The solid line is the inferred local value of the volumetric 
dissipation due to the turbulent cascade. The thick line of triangles represents a term in the 
equation for kinetic energy density (Eq. [1]). Divergence of Fk (left) and divergence of Fp (right) 
are shown. The terms for divergence of flux are important at the boundaries; they smooth the 
distribution of kinetic energy, causing the turbulent velocity to approach the form required for 
Kolmogorov damping. 




6,87 

Fig. 4. — Comparison of terms in the kinetic energy equation, averaged over four transit times 
and over angle, as a function of radius. The solid line is the inferred local value of the volumetric 
dissipation due to the turbulent cascade. The thick line of triangles represents a term in the 
equation for kinetic energy density (Eq. [1]). The terms p'V -u' (left) and dEx/dt (right) are shown. 
The effects of p'V • u' are small and most noticeable at the lower boundary. The change in the 
average kinetic energy is smaller still. 



-17- 



ation vector, and is necessarily large-scale and anisotropic. The large-scale order of this motion is 
destroyed by turbulence, which spreads through all accessible regions. Thus, over all the convec- 
tive region, except right at the boundaries, the flow is made more isotropic, and the Kolmogorov 
damping becomes better approximated by the local isotropic expression 

£ = p{u')I^JId^ (15) 

(per unit volume). After carefully distinguishing between global Kolmogorov damping (Eq. [5|) and 
the local version (Eq. [T5|) , we find that turbulence tends to drive fluxes in such a way as to make 
both valid. 



2.6. The Phase Shift between Damping and Driving 

Having explored the spatial dependence of the KE equation, we now examine its time de- 
pendence. The flow, while wildly fluctuating, has an orderly statistical behavior. To see this, we 
examine the behavior of kinetic energy, integrated over the volume of the convection zone, as a 
function of time. The integral buoyancy flux, qint = /(g • ^ ^)dr, has dimensions of velocity cubed. 

It is convenient to plot kinetic energy and (^j^ on the same graph for detailed comparison. Figure [5] 
shows the time behavior of damping and driving terms. There are two types of time dependence: 
(1) a set of bursts (pulses), and (2) a secular increase in amplitude (related t o the lack of total 



balance between nuclear heating and neutrino cooling over the convective zone (jArnettI (|l996l ). see 
Ch. 10). 

For a static steady state, these two curves would be almost identical. If for simplicity we assume 
a planar convective region and neglect the smaller terms in the KE equation (Eq. [1]) in favor of 
buoyancy and damping, we find the simple result Vrms" I^D ~ f Qdr, or £d/^cz ~ Vrms^ / f Qdr. 

The power spectra of both curves have a peak at 89 seconds. Both exhibit strong fluctuations; 
the buoyancy flux precedes the kinetic energy by roughly 20 seconds. We identif'^ this lag with the 
time it takes for the turbulent cascade to react to changes in the flow. This is suggestively close to 
our previous estimate of the turbulent decay time of about 14 seconds. 

Buoyancy is primarily a property of the largest scales, while damping is a property of the 
smallest. The large separation of scales means that they are not tightly coupled (except on average), 
a characteristic of turbulent systems. 

The buoyancy flux reaches a high value before turbulence can stop its rise, so that it overshoots 
the steady state condition. This leads to excessive velocities, which then cascade to excessive 
damping. In our simulations, the lag in dissipation behind buoyant driving aids fluctuations about 
the nominal steady state condition. 



^In implicit large eddy simulations (ILES) like this, the numerical cascade only goes down to the grid scale. 
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Fig. 5. — Turbulent velocity Vrms'^ = "^Eturb/^cz (solid) and q^^j^ (dotted) in the convection zone, 
versus time. Integral buoyancy flux is qint = Jqz ^^"""^ ^^"^ units of velocity cubed. Power 
spectra for both variables peak at 89 seconds; a transit time is 51 s. The kinetic energy lags the 
buoyant flux by roughly 20 seconds. The average turbulent kinetic energy in the convection zone 
is about 0.64 x 10^^ ergs, with significant fluctuations about that value. 
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This "boom and bust" cycle is reminiscent of prey-predator relations and the logistic map (|Mav 



19761 ). which also have chaotic behavior. We find that an iterated map of the delay in turbulent 
damping is quantitatively inadequate to drive the fluctuations. They are driven by nuclear burning, 
as our simulation (described below) of turbulent decay implies. It will be interesting to explore 
whether this fluctuating behavior depends upon resolution (we expect it to be dominated by the 
largest scale eddies, and weakly dependent on resolution). 



Consider the average conditions from 200 to 800 seconds in the simulations. The average level 
rrbu 
we find 



2 /3 

of turbulent kinetic energy is about 0.7 of q-^^^ . If we ignore the smaller terms in the KE equation. 



aD = io/icz = 1.54 y" Qdr ~ 0.89. (16) 

This is roughly what we get from Table [H so we conclude that to our accuracy, the dissipation 
length in our simulations is roughly the depth of the convective zone, consistent with Kolmogorov 
theory. The appearance of icz here is due to integration over the convection zone, not to any 
assumption about the nature of the largest eddies. 



2.7. The Decay of Turbulence 



The decay of turbulence is a complex topic (lBatcheloiill960l ). Here we will utilize it to provide 



independent estimates of the damping length. These estimates will differ in detail from those of 
driven convection because the flow details change, and the dissipation is a function of the flow 
properties. Nevertheless the sizes of the damping lengths are found to be comparable. 

In the left panel in Fig. [6l we show an independent simulation: the rate of decay of turbulent 
kinetic energy, after oxygen burning and neutrino cooling are artificially turned off at 439 seconds 
(in the middle of the previous simulation). 

If we neglect the small terms for energy flux escaping the convective zone, we may express 
the kinetic energy equation (Eq. [T|) more concisely as D = dK/dt — B, where D is the inferred 
dissipation. After 30 seconds, the buoyancy driving term B becomes small, and dK/dt tracks the 
turbulent dissipation D. Notice that dK/dt is slightly below D, possibly because of entrainment 
of stable matter at the convective boundary which requires energy. If we globally fit the inferred 
damping term by ek = —Mcz Vrms'^/^D, we have in = 4.16 x 10^ cm = 0.97 Icz- 

The right panel in Fig. [6] shows a fit to the kinetic energy, starting 40 seconds after the 
burning was switched off, and the fossil buoyancy had died away. The kinetic energy may be 
represented analytically if dK/dt = D, so that this late part of the damping curve is reasonably 
well approximated by a lower damping length, ijy = 2.6 x 10^ cm = 0.61 icz- 

This independent simulation supports the identification of sub-grid damping with isotropic 
Kolmogorov damping, and the representation of damping by the global expression (Eq. [5|), and 
in ~ (0.6 to 1.0)£cz. 



20 40 60 80 20 40 60 80 

seconds seconds 

Fig. 6. — (left panel) The rate of decay of turbulent kinetic energy after oxygen burning and 
neutrino cooling are turned off. The dominant terms are shown: B = J WsdV is the increase 
of kinetic energy due to buoyant driving, dK/dt is the time derivative of the total kinetic energy 
in the convection zone, and D = J exdV is the inferred dissipation from the turbulent spectrum 
{dK/dt — B). We assume that little kinetic energy escapes the convection zone. The solid line 
shows the damping rate approximated by ek = —MczVrms^/^D with Id = 4.35 x 10^ cm. After 
30 seconds, the buoyancy has died away but damping remains, (right panel) Kinetic energy as a 
function of time (solid line), and fits of the damping after 40 seconds for ^/^/lO^ cm = 2.1, 2.6 
(best) and 3.1. In this phase there is no driving to create anisotropy, and the flow becomes more 
isotropic on average. 
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Table 3. Decay Lengths for Different Flows 



Flow 


OiD = t-D/icz 


Quasi-steady Oxygen Burning 


0.89 


Fossil buoyancy in decay 


0.97 


Pure decay 


0.61 
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These and the previous estimates for ij:, seem to vary beyond the statistical accuracy of 
the analysis. Replacing the complexity of turbulent dissipation by the simple expression ek = 
—MczVrms^/^D-, with constant is a vast change. This quantity Id is a property of the flow, 
and indirectly of the problem b eing addressed. In particular, the pure decay is more isotropic than 
the driven cases (lOgilvid l2003l ). As Table [3] indicates, the ratio (-d/^cz is different for different 
flows, but of order unity. Further invesitgation with different simulations and different physical 
parameters is needed to understand more precisely the general behavior oIId- 



3. General Equations for Stellar Evolution 

How is this approach related to the standard equations of stellar evolution? 



3.1. Internal Energy Equation 



In MLT a connection is made between superadiabatic gradient AV and convective velocity. 
This a llows the convective enthalpy fl ux to be expressed in terms of AV (see lKippenhahn &: Weigert 
(11990), § 7. lciavtonl(|l98,^ ). § 3-5, and lHansen fc Kawaleil ([l99J), § 5). The turbulent kinetic energy 
equation (Eq. [T]) can perform the same function, if the convective velocity is identified with the rms 
velocity {u'"^)^ at a given radius, and used in the enthalpy flux in the same way. This replaces the 
parameterization used in MLT with a physical constraint. 



We need to rewrite the total energy equation (A. 6 in iMeakin &: ArnettI (|2007bl )) in the form 
of the internal energy equation used in stellar evolution. To get the internal energy equation, we 
subtract the kinetic energy equation (A. 12) from the total energy equation (A. 6), and find 



dt{pEi)+V-{MpEi+Pa)) 
-V- 



Fi + F, 



(p'V • u') 



+eK + (pouo • g) + (pe). 



(17) 



The le ft-hand side is just p{dEi / dt+p^dV / dt) in lagrangian (co-moving) coordinates (jLandau &: Lifshitz 
(|l959l ). see Ch. 1), defining po = 1/^, so that 



d{Ei)/dt + {podV/dt) 
1 



-V- 



Fe + Fr 



-{^V ■u') + sk/po 
po 

+(uo-(g- — Vpo)) 
Po 



(18) 
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The first two lines are the usual terms. Here Fr is the heat flux due to radiative diffusion, and in 
this frame the flux of inter nal energy Fi becomes the enthalpy flux carried by convection Fg (see 
Tennekes &: LumlevI (|l972l ). p. 33, and our Eq. [T2|) . The nuclear heating term e includes neutrino 
cooling. 



We w ill rewrite Eq.fTSlin more familiar notation (see lArnetti (| 19961 ) , IClaytonI (119831 ) , iHansen &: Kawalei 



UBM), or lKippenhahn k Weieertl (|l99m l. 

dE/dt + PdV/dt 



1 
P 



V • [Fe + Fr] + e 



-(— V-u') +e/</po 
Po 

+(uo-(g- — Vpo)) 
Po 



(19) 



The flrst line contains the usual formulation. The new terms are: (^V • u') which represents the 
compressional work done by pressure fluctuations (which also appears in the kinetic energy equation 
and we have seen to be small here), £k/pq which is the deposition of heat by the Kolmogorov 
turbulent cascade, and (uo-(g — ^Vpo)); which is zero in hydrostatic equili brium without rotation 
or ex pansion, and drives meridional circulation for rotating, radiative stars (|Tassoullll978l : IClavton 
1983). 



The ek I Po term allows turbulent kinetic energy to do dissipative heating, a new effect not in 
conventional stellar evolution, and it is not guaranteed to be negligible. Through Eq. [T] it couples 
the divergence of the kinetic energy flux (including rotational shear) and the wave energy flux to 
the internal energy. As an internal energy source term, it can generate entropy and cause mixing 
even in radiative regions. For the Sun, this term would give rise to heating in the photosphere, 
chromosphere and corona by shocks and wave motion (pressure and Alfven waves), for example. 

Equations [T] and 1191 represent an extension of turbulent convection theory for stellar evolution, 
in which (1) the algebraic rela tion for conve ctive velocity and superadiabatic gradient is replaced 
by a differential equation (see ISpiegell (jl972l )). and (2) the Kolmogorov cascade is explicit in the 
formulation. Notice that both space and time derivatives appear, making the system nonlocal 
and time dependent, unlike MLT. These derivatives allow the treatment of boundary dynamics in a 
physical way. We will explore speciflc implementations in a subsequent paper, including entrainment 
at convective boundaries. 

Let us now consider some simpliflcations in order to clarify the meaning of these equations. 
In particular, we assume time invariance (so the time derivatives are small on average), no overall 
background motion (uq = 0), and little work done by pressure perturbations on the velocity 
perturbations. As shown previously, these are reasonable approximations except in extreme cases. 
Then Eq. [1] becomes 

(p'g • u') - V • (Fk + Fp) = EK, (20) 
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and Eq. [19] becomes 

pe - V • (Fe + Fr) = -EK (21) 

The coupling of internal energy and turbulent kinetic energy occurs through the Kolmogorov cas- 
cade, which creates internal energy by damping kinetic energy. Eliminating £k, we have 

pe = -{p's • u') + V • (Fe + Fp + Fr + Fk). (22) 

In a steady state, the nuclear energy generation must balance the divergence of fluxes out of the 
region and supply the work needed to maintain the convective flow. This is different from the usual 
formulation used in stellar evolution in that there are new terms {—{p'g • u') + V • (Fk + Fp)), 
which combine to equal ek, the dissipation rate due to the Kolmogorov turbulent cascade. If the 
turbulent velocity is nonzero, these corrections are also nonzero, leading to the conclusion that 
the standard formulation of stellar evolution is wrong in neglecting turbulent heating. The motion 
implied by convection will carry kinetic energy, drive pressure and gravity waves into radiative 
regions, and give local microscopic heating as it dissipates in a transit time, effects that should no 
longer be ignored in stellar evolution. Through its dependence on for a given convective enthalpy 
flux, the kinetic energy flux is dependent upon the equation of state. 

For a low Mach-number flow and radial coordinates, 

(p'g • u') = PrgFe/CpT, (23) 
so that V • F — > dL/dm, then Eq. [22] becomes 

e = - — Le + ^{Le + Lp + Lr + LK), (24) 
nip dm 

where Vq = PrPV/CpT, nip = ^Trr^pHp, and Hp is the pressure scale height. For an ideal gas, 
Vq — > (7 — l)/7, or 0.4 for 7 = 5/3. For gases with a specific heat at constant pressure which 
is large compared to the specific heat at constant volume (such as partially ionized plasma or 
electron-positron plasma), Va is smaller (see below). 



4. Comparison to Other Simulations 
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Table 4. A Comparison of Parameters from some 3D simulations 



Reference 






aT 






EOS 


Bnd. 


Zones 


Meakin k Arnett f 2007b) 


2.0 


0.73 


0.70 


1.22 


e-pair 


yes 


4.0(6) 


Porter & Woodward 


(2000) 


4.5 


2.04 


0.80 


2.70 


ideal 


grid 


6.7(7) 


Chan & Sofia 


(1989) 


4.8 


1.05 


0.83 


2.16 


ideal 


grid 


3.6(4) 


Kim at al. 


C1995) 


6.0 


1.42 


0.85 


2.16 


ionize 


yes 


3.3(4) 


Chan & Sofia 


fl996) 


6.8 


1.30 


0.68 


2.60 


ideal 


yes 


4.8(5) 




MLT 




a/2 


1.0 


a/V2 
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eral groups have compared 3D simulations to MLT predictions ( 


Chan k Sofialll989 


1996: 


Kim et al. 


1995 


1996: 


Meakin & ArnettlboOTb : 


Porter & WoodwardlboOO : 


Porter. Woodward, k JacobslboOC) 


and published sufficient detail to allow easy quantitative comparison. All agree that 
what successful, but derive different values for some of the MLT parameters. Thit 
these parameters may not be universal, but a function of the conditions of the simu 
that our theoretical analysis may be able to put them on a common basis. See also 


MLT is some- 
i suggests that 
ated flow, and 
Abbett. et al. 



(l20o3)^ who have 



(|l997l ): iLudwig et al. I (jl999l ): iTrampedach. et al. I (119991 ): Brandenburg, et al. 
also compared simulations to MLT. 

These simulations are not a homogeneous set, so th at global comparisons on this data set must 



be taken with caution. iPorter &: W oodward hoo d) andlMeakin fc Arnettl (I2nn7bl') used PPM codes 
while the Yale group (|chan fc Sofia . 1989. . .1996 : iKim et al.l 11991 " iRobinson et al.ll2004l'l used com- 



pressi ble vi scous codes with sub-grid scale modelling (and much lower resolution). IChan &: Sofia 
(|l996l ) and iMeakin fc Arnettl (|2007bl ) had the convection zone bounded by st able regions while 
the others d id not; these boundary conditions give rise to new phenomena (jMeakin fc Arnett 
20061 . l2007al ). The deeper convec tion zones develo p ed ni or e anisotropic fiows, and veloc i ties a p- 



proaching the sound speed (e.g., ICattaneo. et al. I (|l99ll ): IWoodward. Porter, fc Jacobsl (j2003l )). 
Porter fc Woodward! (|2000l ) carefully attempted to compensate for these effects, the Yale group 
use d damping to tarn e them, and the soft equation of state and shallo wer depth made them small 
for iMeakin fc Arnettl (j2007bl ) (see below). IPorter fc Woodward (j2000l ) found that they needed to 
shift the a pparent mixing length by about 30 percent, down from a = 3.53 to a = 2.68. The sim- 
ulations of Meakin fc Arnettl ( 2007a ) includ ed an oxygen-bur ning shell (with an e l ectron -positron 
equation of state, see Table [5]), and those of iKim et al.l (jl995l ) and iRobinson et al.l (l2004l ) included 
a solar photosphere (with strong ionization effects in the equation of state). Our simulations are 
driven by nuclear heating at the bottom of the convection zone; the others are drivn by cooling at 
the upper boundary (like the Sun). 

Table [4] summarizes the inferred parameters. The entries are ordered in increasing depth 
of the convective zone as measured in pressure scale height units (icz/Hp), which corresponds to 
increasing asymmetry in the vertical direction. The choice of equation of state (ideal gas, e~e''~-pair 
gas, ionized plasma) and of the boundary conditions affect the simulations. Even the definition of 
the depth of the convective zone may be modified depending on whether the grid includes the stable 
bounding region ("yes") or n ot ("grid"). The bottorn line in the table gives the traditional MLT 
values for several parameters (jKippenhahn fc Weigertlll990l ). The last column gives the number of 
zones on the computational grid. 
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4.1. Convection Parameters 



We no w examine each o f a set o f important convection parameters (see iPorter Sz Woodward 



([20001) and lMeakin &: ArnettI (j2007bl ) for details). These parameters reflect the various uses of the 
MLT parameter a, and are a convenient and concise way to compare the simulations of compressible 
turbulent convection. 



In 3D simulations, a correlation is found between rms temperature fluctuation {T')rms = {T'"^) 2 
and the superadiabat ic gradient AV = V — Vg — Vx , using the conventional notation in astro physics 
of V = d\nTld\nP dKippenhahn k WeigertlEoOol : iHansen k Kawalei]ll994l : lciavtonlll983l ). Here, 
Vq = (91nT/(91nP, taken at constant entropy and composition, and Vx = is the remaining 
part due to possible compositional change. Then, 



{T')rms/To = arAV, 



(25) 



on average, over time and over the volume of the convection zone. Consider low Mach-number flow 
so that p'/po = -Pt [T'/Tq). 



In a single convective roll (e.g., iLorenj (jl963l ): iTrittonI (|l988l )). T' is the temperature per- 
turbation amplitude at the horizontal mid-plane. In the vertical mid-plane of the roll, T' = 
[{dT/dr) — {dT / dr)ad]^cz is the corresponding amplitude. Assuming the amplitudes are compa- 
rable, 

T'/To = {ecz/2Hp)AV, (26) 

implying that ax ~ iczf^Hp, and not a constant. In MLT, ax = a/2, which is a particular choice 
for the assumed flow. 

In the simple picture of a single convective roll, {T')rms/TQ can give a buoyancy torqu e, while 
AV cannot, because t he gr avitational acceleration vector is directed radially downward (jTritton 
19881 ) . In the iLorena (jl963l l model of thermal con vection, the difference bet ween the two is inti- 
mately connected to the onset of chaotic behavior. iMeakin fc ArnettI (|2007bl ) find ~ 0.73 (see 
their Fig. 17), so that 

ar = 0.73 icz/2Hp. (27) 

There is considerable variation in the values of ay shown in Table HI with a tendency to increase 
for deeper (more stratified) convection zones. 

Notice that for two convective rolls, one atop the other, we would expect the characteristic roll 
size to change {i « £cz — > ^czf^, so ax — > aT/2, approximately). This explicitly shows how the 
convection parameters can be a function of the properties of the flow itself. 

Figure [71 left panel, shows the behavior of ut as a function of th e depth of the convection 
zone, for the computations listed in Table [H The two PPM calculations ([Meakin &: ArnettI (l2007bl ) 
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and lPorter &: Woodward! (|200d )) agree with the scaling indicated in Eq.[271 The calculations using 
the compressible viscous equations with sub-grid modelling for dissipation all lie below Eq. [271 In 
order to give some idea of the depth needed in simulations of stellar convection zones, the x-axis 
in Fig. [7] is marked from zero to the depth of the solar convection zone (20 pressure scale heights). 
None of these simulations describe such an extremely anisotropic case. 



4-1-2. oe and ax 

The enthalpy flux is 



F, = poCpT'u' 

(XEPoCp{T')rms {u')rms, (28) 



where iMeakin &: ArnettI (|2007bl ) find aE = 0.70 it 0.03; in Tabled! a^; is relatively constant among 
the simulations (the total range is about 20 percent). It is not ruled out that oe might be a 
universal constant, or at least slowly varying. Note that q^; is just the correlation coefficient for T ' 
and u'. It seems that a_E is not sensitive to the Prandtl number Pr. iPorter &: Woodward! ()2000l ) 
have a different Pr (ours is Pr ~ 1) and get an aE similar to ours. 

Using Eq. [25! this becomes 

Fc = a£aT/OoCpro(n')rmsAV. (29) 
Similarly, if the kinetic energy flux is 

Fk = \{pu'^0, (30) 

we may define 

aK = {pu'^<)/{p){Kffl\ (31) 
so that Fk = ^(p) ((u^)^)^/^. Note that the sign of ax can be negative. 



4.1.3. a„ 

In lMeakin &: ArnettI (l2007bl ) . the correlation between convective velocity (squared) and super- 
adiabatic excess AV is written as 



(n')Ls = (a^/4)5/3T^pAV. 



(32) 



In MLT, a similar expression is defined, with a^/8 replacing If we have local balance between 

buoyancy driving and turbulent damping, 



{p'gu') « {p){u')l^Jl 



D- 



(33) 
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using Eq. [131 Eq. [25l and Eq. [28l we have 




rms 



{lDaT0lE)l3Tg^^ ■ 



(34) 



Comparing this to Eq. [32] we have 



all A = aTOtE^D/Hp. 



(35) 



Using Eq.[27|and a^; = 0.70, 



= 1.22 (^cz/2//p)(^d/0.9^cz)5 



(36) 



In Table [H we see that a„ is variable, and tends to increase with increasing depth of the convective 
zone. This is shown explicitly in the right panel of Fig. [71 The two PPM calculations are in 
excellent agreement with Eq. [36l taking ip) = 0.9 icz {cud = 0.9), while the three compressible 
viscous calculations lie below it. 

Using Eq. [29l we have 



Fe = [aEaTa„/2)poCpn^/gPTHp{AVf/\ 



(37) 



Using Eq,[271 Eq,[36l = OMcz and ue = 0.70, we have 



OEaraJl = 0M2{ecz /Hpf , 



(38) 



for the factor in Eq. [37l 



4.2. The electron-positron Plasma 
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Meakin-Amett(2o67) 
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Fig. 7. — Predicted 2ar = aA,T (top panel) and V^a^ = a a, 7; (bottom panel), a s a function of 
depth of the convective zone in units of Hp. This scaling (jMeakin &: Arnettll2007bl ) gives "alphas" 
which are comparable to the MLT values and each other. Some of the error bars are large; new 
simulations are needed to determine how well such results follow a single curve (Meakin &; Arnett, 
in preparation). In the limit of small depth, the mixing length m ust be no larger th a n the d epth 
itself; hence the point at zero depth is an analytic result. The iMeakin &: ArnettI (|2007bl ) and 
Porter &: Woodward! ((20001) points and the origin are close to colinear. It appears that both ot 
and are functions of depth of the convective region, and not universal constants. 
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Table 5. Thermodynamic parameters for Entropy S/TZ = 4.623 





P6 


In P-max / P 


Cp/n 


/3t 


PV/CpT 




2.51 


1.580 


0.0 


15.36 


3.912 


0.0609 


0.2382 


2.31 


1.225 


0.348 


14.53 


3.757 


0.0638 


0.2397 


2.11 


0.930 


0.724 


13.62 


3.592 


0.0673 


0.2417 


1.91 


0.690 


1.134 


12.67 


3.421 


0.0714 


0.2443 


1.71 


0.498 


1.583 


11.73 


3.257 


0.0762 


0.2483 


1.51 


0.348 


2.080 


10.83 


3.106 


0.0814 


0.2529 
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The thermodynamic character of the electron-positron plasma in the oxygen burning shell is 
significantly different from an ideal gas, an effect which must be taken into account in comparisons 
to simulations which use differe nt equations of sta t e. Th is is illustrated in Table [5l we use the 
Helmholtz equation of state of iTimmes fc: Swestvl ([20001) . The entropy in the oxygen burning 
convection zone is S/IZ ~ 4.6 in dimensionless units, where TZ is the gas constant. The zone is 
about 2 pressure scale heights deep (see column 3). For an ideal gas, the specific heat at constant 
pressure is 2.5TZ; Cp/lZ is much larger for the plasma, ranging between 10 and 16 (column 4). 
Column 5 gives the value of /3t, which is unity for an ideal gas, but lies between 3 and 4 for the 
plasma. Column 6 gives the ratio of PV/CpT = pq/ pqCpTq, which is 0.4 for an ideal gas. The 
ratio of the buoyancy flux to the enthalpy flux (Eq. [T3]) is proportional to Vq = I3tPV/CpT. The 
same factor appears in the ratio of the kinetic energy flux to the enthalpy flux. For the ideal gas 
Va = 0.4, but is smaller for the plasma (column 7). For the same convective enthalpy flux, the 
pair-plasma has a lower kinetic energy flux, so the velocity scale is lower. 



4.3. Direction of the Kinetic Energy Flux 

The kinetic energy flux is averaged over angle at a given radius, and averaged over two convec- 
tive turnover times. As Fig. [5] shows, this time span covers significant dynamic behavior. Roughly 
speaking, an unstable configuration forms, becomes dynamic, reforms, becomes dynamic again, and 
so on. Over this turnover timescale the nuclear burning provides the energy necessary to drive the 
turbulent kinetic energy (er ^(u'^)). 



The convective instability is closely related to the Rayleigh- Taylor instability (IChandrasekhar. S 



1961), which has been produced dramatically in high energy-density plasma experiments (iRemington et al 
19991 ). i.e., under star-like conditions and well into the nonlinear growth regime. The characteristic 
behavior is the rise of mushroom-shaped higher entropy plasma and the fall of spike-shaped lower 
entropy plasma. If we consider a closed volume, these motions are accompanied by slower return 
motions which maintain mass conservation. The kinetic energy flux scales with velocity cubed, and 
so is dominated by the fast moving mushrooms and spikes. In a symmetrical system we will have 
an upward kinetic energy flux (from the mushrooms) in the region above the horizontal midplane 
of the volume, and a downward kinetic energy flux in the region below. If we average over several 
cycles (turnover times) the kinetic energy flux with be dominated by these motions, being positive 
(upward) above the midplane and negative below. This is qualitatively similar to the simulation 



results of lMeakin ArnettI (j2007bl ) (see their Fig. 21 and Fig. 22). 



This simple picture is complicated by an up-down asymmetry due to stratification (hydrostatic 
structure). The depth of the convection zone in pressure scale heights, (.cz/Hp = ln(Pf,of /Ptop), 
is a direct measure of the up-down asymmetry. Upward flows move into regions of lower pressure, 
and expand; downward flows move into regions of higher pressure, and are compressed. The 
smaller area of the downflows implies higher velocities relative to a coordinate frame containing 
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the convection zone (a lagrangian frame). This will favor dow nward (negative) kinetic energy 
fluxes. The neglect of kinetic energy fluxes (|Bohm-Vitensdll992l ) corresponds to the limiting case 
of a shallow conve ction zone. For deep convection zones, there is a strong bias in favor of fast 
downward plumes (jNordlund k. Steinll2000l: IStein &: Nordlundlll998l '). and these dominate the flow 
for simulations with icz/Hp > 4 or so. In lMeakin fc ArnettI (|2007bl ). which has a depth Icz I Hp = 
2, the convective zone is skewed, so that the surface, which separates the posi tive and the negative 
kineti c energy fluxes, moves nearer to the bottom of the convection zone. For lPorter Woodward 
([20001), where icz/Hp is larger, the positive kinetic energy flux is overwhelmed, and the direction 
of the kinetic energy flux is opposite to that of the much larger enthalpy flux. We expect this to 
be a general property of deep convection zones. 

There is a natural limit to the depth of convection zones expecte d in active burning regions. 
Entropies in active burning regions vary relatively slowly (jArnettI (119961 ). Ch. 10), so that the depth 
of a convection zone implies a value of the temperature ratio between bottom and top. In Table O 
the convection zone extends down almost to neon burning temperatures (T ~ 1.5 x 10^ K). Further 
extension will entrain new fuel into the oxygen convective shell, which will burn at the top of the 
convection zone, choking the flow. For the last stages (C, Ne, O and Si burning), the burning 
temperatures for different fuels are fairly close together, implying that their related convection 
zones will tend to be relati vely shallow. They are a lso close together in radius, raising the issue of 
interactions between them (jMeakin &: ArnettI l2006l ). 



Here icz/Hp = In Phot /Ptop = Tpj^^Thot/Ttop- For example, for advanced burning stages 
or radiation dominated regions, 7 ~ 4/3, and 7/(7 — 1) = 4. For helium and hydrogen burning, 
T{He) /T{H) ~ 13 so icz/Hp < 41nl3 ~ 10. A helium burning convective zone will not be deeper 
than about 10 scale heights unless the overlying matter is devoid of hydrogen. For carbon burning 
the corresponding depth is icz/Hp < 3, unless devoid of H and He fuel. Surface convection zones 
may extend down to the hydrogen burning regions. Very roughly, lnT{H)/T(, 8, so icz/Hp < 32, 
using the structure of an n = 3 polytrope. 



4.4. Magnitude of the Kinetic Energy Flux 
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Table 6. Comparison of some 3D simulations 



Reference 








Meakin & Arnett f 2007a) 


2 


0.24 


-0.0f8 








+0.014 


Porter & Woodward (2000) 


4.5 


0.40 


-0.3 


Cattaneo. at al. (199fl 


5 


0.40 


-0.35 


Chan & Sofia (1989) 


4.8 


0.40 


-0.35 


Chan & Sofia ff996) 


6.8 


0.40 


-0.50 
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Table [6] gives the ratio of kinetic energy flux to enthalpy flux for several 3D simulations. This 
ratio is much smaller in our simulations (by a factor of 35 to 50) than in the others, and Fk changes 
sign in our convection zone. 



As in (iMeakin &: Arnettll2007bl ) we can write the KE to enthalpy flux ratio according to mixing 



length relationships, 

Fk pvl/2vc ( al 



Va (39) 



Fc pCpT' Vc \8aT, 
and then balance between buoyancy driving and turbulent damping through Eq. [35] gives 



Fk aooiE ( Icz 



Fc 2 \Hp 



Va. (40) 



See Eq. 3.14 of iPorter &: Woodward! (j200d ). which uses a gamma-law equation of state to 



generate the sum of kinetic and enthalpy fluxes, and implies an equivalent flux ratio. 

The ratio of for the ideal gas to that of the electron-positron gas gives a factor of 1.6 
or so. The ratios of the depth of the convection zones give another factor of ~ 5. While these 
considerations illustrate the role played by the depth of the convection zone and go some way 
towards explaining the trends in KE flux it is also important to consider the geometry of the 
driving region, which relates to how well a local balance between buoyancy driving and turbulent 
damping is achieved as assumed in Eq. [35l In particular, the length scale la over which buoyancy is 
imparted to the stellar plasma through either heating at the base of the convection zone or cooling 
at the top can affect the KE to enthalpy flux ratio dramatically. This may be understood simply: 
the flow depends both upon the geometry of the convective domain and upon the way in which the 
fluid is stirred. At present we have at least two basic patterns, a mostly negative kinetic energy 
flux for deep convective zones driven by surface cooling (most stellar surface convection zones) and 
a more complex positive-negative convective flux for shallower zones driven by nuclear burning. 
More simulations are underway to clarify this issue (Meakin and Arnett, in prep.). 



4.5. Saturation of Kinetic Energy Flux 



Are there limits to the linear rise in energy of convection that is implied by Fig. [7p Deep 
convection zones {icz ^ ^Hp) have strong negative kinetic energy fluxes. For very deep zones, 
the extrapolated kinetic fluxes imply supersonic velocities. Such large velocities would generate 
shock waves, which would increase dissipation, converting kinetic energy into internal energy. The 
rate scales as the velocity difference cubed, which is the same scaling as turbulent dissipation in 
the Kolmogorov cascade (|Bethdll942l : iBorid 120071 ) . This suggests that even if deeper convection 
zones did tend to have increasingly strong velocities, shock dissipation will become significant, and 
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resist the increase with increasing £cz- In this sense, the increase in damping length must 

"saturate" with increasing depths. 

Increased damping may come sooner from another effect: a change in the nature of the eddies 
and the size of the damping length. Physically this would occur as follows: as the deep, fast, 
narrow downflows drive into the convection zone they will give rise to shear instabilities at their 
interfaces, which will lead to mixing with the ambient fluid, and exchange momentum through a 
turbulent viscosity, and eventually completely dissolve into the background. This would give shorter 
dissipation lengths, and would begin to occur before the mach numbers become large enough for 
significant shock formation. 

We expect the PPM calculations for deeper convective zones to "flatten" in Figure [7] as the 
Yale simulations do, but at higher amplitude, due to increased damping with increased depth. This 
hypothesis needs to be tested numerically, which will be challenging. Shock waves may cause errors 
at grid boundaries, deep convection zones will have longer thermal relaxation times, and maintain- 
ing sufficiently wide aspect angle implies many computational zones for adequate resolution, for 
example. 

The simulations of the Yale group (Ichan fc Sofialll98d.[ 



1996 



20041 ) ■ which use fewer zones and a str ong damping (ISmagqrinskv 
dissipation than the PPM simulations ( Porter &: Woodward ]20od : ' 



Kim et al.lll995l . 119961 : iRobinson et al 



19631). appear to have a larger 



Meakin fc ArnettJl2007bl ). This 



is qualitatively equivalent to the saturation discussed above, and may be tested if calculations using 
the Yale code (or equivalent) are repeated at higher resolution and/or lower dissipation. 



5. Some Implications 
5.1. Waves 

The energy flux terms include both a dvective transport a nd wa v es. Here we will recall some 
properties of waves and their generation ( Landau Sz Lifshitz (1959), Press! (jl98ll )). The charac- 
teristic frequencies of convective motion are centered about a frequency tocz ~ Vrms/^cz- The 
convective mach number, {u')rms/C = p' / Po = p' /poC^, is a measure of compressibility of the flow; 
here C is the sound speed. If the interface between convective and radiative zones moves with the 
matter (is volume conserving , on average), it generate s acoustic waves by dipole emission at a lu- 
minosity Lp_mode '^%z (see lLandau &: Lifshita (|l959l ). § 73), or as the Mach number to the sixth 
power. For more vigorous motion, the perturbation may give volume changes, so that the acoustic 
wave generation by monopole emissivity, or Mach number to the fourth power. For subsonic flows 
this channel is closed to significant energy flow, but open as {u')rms approaches the local sound 
speed, as it does in the surface convection zones of many stars, including the Sun, or as it may in 
the stage prior to core collapse in massive stars. 



While the exact power dependence may depend upon specific geometries and degree of interfer- 
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ence, a general result seems to be: the gravity wave channel d ominates over the acoustic channel for 
low mach number flows (as we observe in our simulations, see lMeakin &: ArnettI (|2007al ) for a more 
detailed discussion of both the g-mode and p-mode behavior, including mixed p- and g- modes). 
Both types of waves are generated by convection interacting with stably stratified bounding re- 
gions, and the luminosity of each depends upon both the convective vigor and the impedence at 



the boundary. Such waves ca n propagate into stably stratified regions (jPress 



1981 



1981 



Press fc Rybicki 



Young &: Arnettil2005l ). Generation of gravity waves by conv ective turbulence has become 



an issue in questions of rnixing and angular momentum transport (IGarcia Lopez &: Spruit 



Charbonnel k Talonlll999l : iTalon k Charbonnelll2004l : lYoung ArnettI l2005l ). 



1991 



The establishment of a robust estimate of the turbulent velocity field should improve estimates 
of wave generation, which has implications for mass loss, mixing in radiative regions, coronal heating 
and helioseismology. In particular, the pressure correlation flux takes over the energy carried by 
the kinetic energy flux as the convective boundaries are approached . This gives a direct connection 
between the scale of turbulent velocity and g-mode wave emission ()Presslll98ll ). 



5.2. Rotation and Magnetic Fields 



A closely similar set of mean-field equations are used in the theory of the magnetic reson ance 



instability (MRI) in accretion disks (jBalbus &: Hawlevi Il998l : iPessah. Chan. &: Psaltid 120061 ) . If 
we had included magnetic fields in the MHD approximation, and assumed strong rotation, our 
procedur e would have produced an equation for mechanical energy (our Eg. [1] corresp onds to 
Eq. 17 of balbus fc Hawlevi (jl998l )). and total energy (Eq. A6 of lMeakin fc ArnettI (l2007bl ) to their 
Eq. 27). Projection onto a cylindrical coordinate system, with the rotational axis oriented parallel 
to the total angular momentum vector, would give an angular momentum equation (their Eq. 29.). 
Auxiliary equations provide for magnetic field dynamics, radiation transport, and nuclear burning. 
This underlying similarity provides a way to write a more general set of mean-field equations, of 
which both stars and accretion disks are limiting cases. In turn this allows a systematic evaluation 
of the relative importance of different effects (rotational mixing and convective mixing, for example) 
which are now considered piecemeal. 

Our simulations include the complete set of rotational terms, but the initial conditions im- 
ply that these terms are not exercised except as convectively induced shear. Our simulations 
do not include magnetic fields in the MHD approxim ation, but could be generalized to do so 
( Stone &: Gardineiil2007l : iPessah. Chan. &: Psaltisll2006l ). Because they interact, rotation and mag- 
netic fields should be included together. 



Balbus Hawleyl (|l994l ) have argued that in the stellar case, the weak-field MRI dominates 



over merely hydrodynamic in stabilities, and drives the radiativ e zones (but not convective zones) 
toward solid body rotation. iHeger. Langer. k. Woosleyl (j2000l ) argued the reverse (based on the 
H0iland instability criterion), that convective zones would tend toward rigid body rotation, and 
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radiative zones would tend to have differential rotation. Because the H0iland instability ap- 
plies to neutral fluids, not dense plasma, it is probably not r elevant for stars. Helioseismology 
(jThompson. et al. Ill996l : iHowe. et al. II2005I : iBrandenbur a 120071 ) is showing that while the convec- 
tive zone of the Sun shows differential rotation, the underlyin g radiative zone seems t o be tending 
toward solid body rotation, as the MRI arguments suggest. (jBrowning &: Basil 120071 ) have shown 
that even deep convection zones can generate significant magnetic fields, so that the presence of a 
stable interface is not necessary for field generation. The rotational state of the central regions of 
the Sun probably depends upon the efficacy of a ngular momentum transpor t by processes related 
to fiow of the plasma, including g-mode waves (jCharbonnel &: TalonI Il999l ). as well as magnetic 
fields. 



5.3. Damkohler Number for Burning 



In general, it is appropriate to decompose the equations in temperature (T = Tq + T') and 
composition as well. T he opacity and the nuclear reaction rates are often sensitive to both. 



Meakin fc Arnetti (|2006l ) found fiashing due to oxygen burning in vigorous downdrafts which were 
fuel-rich (the fiashes were too mild to affect the fiow dramatically). None of these effects are in- 
cluded in standard stellar evolution theory. For simplicity we supress this complication for the 
moment; this means that our opacities and reaction rates are to be interpreted as averages over 
fluctuations in these variables as well. Further investigation of this issue, with 3D simulations, is 
desirable. 

For this set of simulations, the heating and cooling times are at least 100 times longer 
than the transit times, so that we are in the regime of small Damkohler number Dm ^ 0.01 
(jZel'dovich. Barenblatt. Librovich. &: Makhviladzd Il985l : I Or an &: BorisI 119871 ) . The release of nu- 
clear energy during a transit time is small relative to the internal energy in the convection zone, but 
comparable to the superadiabatic energy and to the turbulent kinetic energy. The burning drives 
the turbulent motions, but only gives moderate pulses of kinetic energy, as shown in Fig. [5j This 
is unlike the much more complex problem of Type la supernova models, for which Dm is large. In 
our case the necessary averaging over convective cycles does not seem to be a problem. 

This convenient state may not apply to the double shell flash stage for asymptotic giant 
branch stars, in which wave driven mixing and entrainment are likely to co mplicate the issue 
of mixing, and therefore flg ure into the question of s-process nucleosynthesis (IBusso et al.lll999l : 
CampbeU LattanziJbood ). 



6. Summary 



We flnd that our three-dimensional time-dependent simulations of compressible stellar turbu- 
lence (ILES) are well represented by a master equation (Eq. [1]) for kinetic energy which includes 
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dissipation implied by the Kolmogorov cascade. The damping length is found in three independent 

ways, with reasonable consistency, and is the size of the largest eddy, which is approximately the 
geometric linear dimension (depth) of the convcctivc zone. Unlike the mixing length of MLT, it is 
not a free parameter, but a mathematical consequence of the turbulent flow. Balancing turbulent 
buoyant driving with Kolmogorov damping provides a reasonable estimate of the damping length. 
The divergence of kinetic energy and acoustic fluxes is nonzero, and provides the mechanism to 
spread turbulence through the convective zone, and make the Kolmogorov damping a good approx- 
imation. The turbulent flow is highly dissipativc and must be maintained by continual driving; 
this "frictional heating" term is missing from the standard stellar evolution equations, as are the 
kinetic energy and acoustic fluxes. 

Fluctuations in kinetic energy are significant (of order 50 percent), and damping lags driving 
by about a transit time for the convective zone. 

Comparison with some other simulations, which were dramatically different in many respects, 
gives a consistent picture. Turbulent convection is more vigorous for deeper (more stratified) 
convection zones. Turbulent kinetic energies are lower for nonideal equations of state, such as 
partially ionized plasmas and electron-positron plasmas, to the extent that their specific heat at 
constant volume is less that their specific heat at constant pressure. The average flow structure 
changes in a simple way depending upon the depth of the convection zone; deep convection zones 
have downwardly directed flows of kinetic energy, cancelling some of the upward enthalpy flux. 

It appears that extension of this approach, using simulations to define stellar convection algo- 
rithms, can establish a theoretical model of turbulent convection that docs not require astronomical 
calibration, but can be based upon a combination of computer simulations, terrestrial observations, 
and experiments. Efforts to implement this general theory in a stellar evolution code are underway. 

During the course of this project, we lost two friends who were leaders in the field of stellar 
evolution, John Bahcall and Bohdan Paczynski, to whom this paper is dedicated. This work was 
supported in part by NSF Grant 0708871 and NASA Grant NNX08AH19G at the University of 
Arizona, and the ASCII FLASH center at the University of Chicago, One of us (DA) wishes to thank 
the Aspen Center for Physics and the International Center for Relativistic Astrophysics (ICRA) for 
their hospitality, Brian Chaboyer for help with the history of the mixing length implementation, 
Martin Pessah for discussions of MRI physics, Robert Stein for discussion of the effect of the 
continuity equation on flows, Vittorio Canuto for helpful discussions of the philosophy of turbulence 
modelling, Martin Asplund for providing machine-readable copies of solar surface models, and 
Frank Timmes for helpful comments on the draft and for providing access to the Saguaro computer 
cluster. We wish to thank an anonymous referee for insightful comments which improved both our 
presentation and our understanding. 
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A. Decomposition 

We decompose velocity, density, and pressure fields into mean and fluctuating components 
according to 

^p = ^PQ + <^', (Al) 

where = ipQ and {ip') = and the overbar and brackets indicate time and horizontal averaging, 
respectively. For data handling and analysis purposes we consider the fluctuating component of 
the field to be composed of a radial p-mode component and a component due to all other sources 

^' = ^'p + ^'t. (A2) 

This additional decomposition allows us to make a more accurate estimate of the fluctuations 
associated with turbulent convection in the presence of a coherent radial p-mode which is not well 
sampled in the output flies from the simulation. If the radial p-mode contribution were well sampled 
then we would flnd {ip'p) = to the degree that the mode is adiabatic. 

Consider the radial velocity to be composed of the following components 

u = uq + u'p + u'f, (A3) 

where uq is the mean background expansion, u'p is the radial p-mode induced fluctuation, and u'f. 
is the fluctuation due to turbulent convection and other non-radial modal components. Averaging, 
we find 



{u) = {u>p) + {u't)+uo. (A4) 

Because of our poor sampling of the low order radial p-modes which have frequencies comparable 
to the simulation data output rate {5t ~ 0.5 s) we find that the term (u'p) > and contributes a 
significant error to the estimation of u[. In order to correct for this horizontally coherent p-mode 
induced radial displacement, we subtract the horizontally averaged radial velocity component at 
each time step and estimate the turbulence induced fluctuation by 

u'f. = u— {u'p) — no ~ « — {u'p). (A5) 

The latter approximate equality in the above equation is made because of the smallness of the 
background expansion compared to the the r.m.s. velocity fluctuations associated with the turbulent 
convection, uq/u'^ ^ 10~^. 
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The instantaneous fluctuations in pressure and density are calculated according to 



Pt=P - (p) 



(A6) 



and 



Pt = P - (P)- 



(A7) 
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